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The critical behavior of the disordered ferromagnetic Ising model is studied nu- 
merically by the Monte Carlo method in a wide range of variation of concentration 
of nonmagnetic impurity atoms. The temperature dependences of correlation length 
and magnetic susceptibility are determined for samples with various spin concentra- 
tions and various linear sizes. The finite-size scaling technique is used for obtaining 
scaling functions for these quantities, which exhibit a universal behavior in the criti- 
cal region; the critical temperatures and static critical exponents are also determined 
using scaling corrections. On the basis of variation of the scaling functions and val- 
ues of critical exponents upon a change in the concentration, the conclusion is drawn 
concerning the existence of two universal classes of the critical behavior of the di- 
luted Ising model with different characteristics for weakly and strongly disordered 
systems. 
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I. INTRODUCTION 

Analysis of the critical behavior of disordered systems with quenched structural defects 
is of considerable theoretical and experimental interest. Most real solids contain quenched 
structural defects, whose presence can affect the characteristics of the system and may 
strongly modify the behavior of the systems during phase transitions. This leads to new 
complex phenomena in structurally disordered systems, which are associated with the effects 
of an anomalously strong interaction of fluctuations of a number of thermodynamic quanti- 
ties, when any perturbation introduced by structural defects (even in small concentration) 
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may strongly change the state of the system. The description of such systems requires the 
development of special analytic and numerical methods. 

The following two questions arise when the effect of structural disorder on second-order 
phase transitions is investigated: (i) do the critical exponents of a homogeneous magnet 
change upon its dilution by nonmagnetic impurity atoms? and (ii) if these exponents change, 
are the new critical exponents universal (i.e., independent of the structural defect concen- 
tration up to the percolation threshold)? The answer to the first question was obtained 
in l|, where it was shown that the critical exponents of systems with quenched structural 
defects change as compared to their homogeneous analogs if the critical exponent of the heat 
capacity of a homogeneous system is positive. This criterion is satisfied only by 3D systems 
whose critical behavior can be described by the Ising model. A large number of publica- 
tions are devoted to the study of the critical behavior of diluted Ising-like magnets by the 
renorm-group methods, the numerical Monte Carlo methods, and experimentally (see review 
{2!). An affirmative answer has been obtained to the question concerning the existence of a 
new universal class of the critical behavior, which is formed by diluted Ising-like magnet. It 
remains unclear, however, whether the asymptotic values of critical exponents are indepen- 
dent of the rate of dilution of the system, how the crossover effects change these values, and 
whether two or more regimes of the critical behavior exist for weakly and strongly disordered 
systems; these questions are the subjects of heated discussions. 

This study is devoted to numerical analysis of the critical behavior of a diluted 3D Ising 
model in a wide range of concentration of quenched point defects. The fundamental impor- 
tance of the results of this study is due to stringent requirements to simulation conditions 
imposed in the course of investigations; the wide range of linear dimensions of lattices 
(L = 20 — 400) analyzed in this work; the chosen temperature range of simulation close to 
the critical temperature with r = (T — T^/T^. = 5 ■ 10^^ — 10~^, which makes it possible 
to single out the asymptotic values of characteristics; high statistics used for averaging of 
thermodynamic and correlation functions over various impurity configurations; the applica- 
tion of finite-size scaling technique jsj] for processing the result of simulation, which makes it 
possible to obtain scaling function for thermodynamic functions apart from their asymptotic 
values; and application of corrections to scaling for determining the asymptotic values of 
critical exponents. 



II. COMPUTER SIMULATION TECHNIQUE AND RESULTS 



We consider a model of a disordered spin system in the form of a cubic lattice with linear 
size L under certain boundary conditions. The microscopic Hamiltonian of the disordered 
Ising model can be written in the form 



where Jij is the short-range exchange interaction between spins fixed at the lattice sites 
and assuming values of ±1. Nonmagnetic impurity atoms form empty sites. In this case, 
occupation numbers pi assume the value or 1 and are described by the distribution function 



with p = 1—c, where c - is the concentration of the impurity atom. The impurity is uniformly 
distributed over the entire system, and its position is fixed in simulation for an individual 
impurity configuration. We consider here disordered systems with spin concentrations p = 
0.95, 0.80, 0.60, and 0.50. 

To suppress the effect of critical slowing down and correlation of various spin configura- 
tions, we used the single-cluster Wolf algorithm, which is most effective in this respect jj, . 
A Monte Carlo step per spin (MCS) was assumed to correspond to 1020 rotations of a Wolf 
cluster depending on the linear size of the lattice being simulated, the spin concentration 
of the system, and the closeness of the temperature to the critical point. The stabilization 
of thermodynamic equilibrium required 10"^ Monte Carlo steps, and 10^ steps were allotted 
to statistical averaging of quantities being simulated for a given impurity configuration. To 
determine the average values of thermodynamic and correlation functions, averaging over 
various impurity configurations was carried out along with statistical averaging (averaging 
was carried out over 3000 samples for p = 0.95, over 5000 samples for p = 0.80, and over 
10000 samples for p = 0.60 and 0.50). 

In the course of simulation of various spin systems, correlation length and susceptibility 
Xl were carried out on lattices with a linear size L in accordance with the following relations: 
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where (xi^j, X2,i, xs^j) are the coordinates of the i-th lattice site; angle brackets indicate 
statistical averaging over Monte Carlo steps, and the bar indicates averaging over impurity 
configurations. The temperature dependences ^l(T) and Xl(T) were determined in the 
temperature interval r = 5 ■ 10~^ — 10~^ for samples with p = 0.95 and a linear size in the 
range of L = 20 — 400. For samples with the remaining spin concentrations, temperatures 
were chosen in the interval of r = 10^^ — 10^^ for values of L ranging from 20 to 300. In 
computer simulation, the value of L^ax for each temperature was limited by the lattice size 
for which the correlation length and susceptibility of the system attained their asymptotic 
values. 

In accordance with the results obtained in 0, 5] and the results of our investigations, 
the chosen simulation conditions ensure equilibrium values for measurable thermodynamic 
quantities for all lattice sizes and spin concentrations studied here since the autocorrelation 
times for magnetization and energy turn out to be not longer than ten Monte Carlo steps 
per spin even for chosen temperatures closest to the critical temperature (with allowance for 
the number of turns of the Wolf cluster taken as a step). 



III. METHOD OF FINITE-SIZE SCALING 



It is known that the second-order phase transition considered here can be manifested only 
in the thermodynamic limit, when the volume of the system and the number of particles in 
it tend to infinity. To determine the asymptotic values of thermodynamic quantities A{T) 
exhibiting an anomalous behavior near the critical temperature from their values Al{T) 
determined on finite lattices, the concepts of the scaling theory concerning the generalized 
uniformity of thermodynamic functions in the critical region relative to scale transformations 
of the system are widely used. These concepts formed the basis of various methods of finite- 
size scaling. Here, we apply the method proposed in and tested by the authors in analysis 
of the results of simulation of the critical behavior of 2D and 3D pure Ising models. 

n 

The idea of this method [3] is that, in accordance with the scaling theory, the size depen- 
dence of a certain thermodynamic quantity defined on a finite lattice in zero magnetic 




FIG. 1: Scaling functions for (a) correlation length and (b) susceptibility obtained at various 
temperatures for a system with p = 0.95 using the approximation polynomial in x. 

field can be presented in the critical region in the form 

Al{t) = L'/''Usl{t)), sl{t) = L/Ur), (5) 

where S is the critical exponent for the thermodynamic quantity A{t) ~ t^^. Taking into 
account the fact that the correlation length in the critical region behaves as ^{t) ~ r^'^, we 
can write 

L'/-^ = A(r).f (r). (6) 
Then expression can be written in the form 

A^ir) = A{t)Fa{sUt)), (7) 

where the relation between scaling functions J'a and is defined in the form of the relation 

FA{sL{r)) = 4\t)Us,{t)). (8) 

If correlation length ^ plays the role of quantity A, Eq. ([7]) defines ^l{t)/L as a function 
of only one variable ^{t)/L. This leads to a relation that makes it possible to find the 
asymptotic value of any thermodynamic quantity in terms of directly measurable values of 
Al and the scaling function of xl{t) = C,l{t)/L, 

A{t) = Al{t)/Qa{xl{t)), (9) 
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FIG. 2: Scaling functions for (a) correlation length and (b) susceptibility obtained at various 
temperatures for a system with p = 0.95 using the approximation polynomial in exp(— 

where function Qa{xl{t)) is defined by the expression 

Qa{xl{t)) = FA{f[\xL{T))). (10) 

ScaHng function Qa{xl), defined in the interval < xl < Xc, where Xc is the value of 
the argument independent of L in the critical region, must satisfy the following asymptotic 
conditions: lira Qa{x) 1 and lim Qa — * 0. 

To satisfy the asymptotic conditions, we chose, analogously to |3|, the scaling function 
for susceptibility and correlation length in the form of a polynomial dependence of x, as well 
as of exp (— : 

QAix) = 1+CiX + C2X^ + Cj,X^ + CiX^, (11) 

Qa{x) = 1 + cie-i/" + C2e-2/- + c3e-3/- + C4e-^/-, (12) 

with coefficients c„ selected for each temperature T using the least squares method. 
Here, we implement the following scheme of finitesize scaling. 

1. For an arbitrary value of tq in the critical temperature range, the values of Al{tq) and 
x(L, To) = C,l{tq)/L are measured for lattices with increasing size L. 

2. The thermodynamic value of quantity A{tq) is determined as the value of Al{tq), 
which is found to be independent of L within the error of measurements. 
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FIG. 3: Scaling functions for (a) correlation length and (b) susceptibility obtained at various 
temperatures for a system with p = 0.50 using an approximation polynomial in x. 
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FIG. 4: Scaling functions for (a) correlation length and (b) susceptibility obtained at various 
temperatures for a system with p = 0.50 using the approximation polynomial in exp(l/2;). 

3. The results of measurements for Al{to)/A{to) are processed by the least squares 
method to determine the corresponding functional form for scaling function 

QAix{L,To). 

4. The procedure is repeated for other values of r in the range of r ~ 10~^ — 10~^. 

5. Averaged scaling function Q^^r is determined on the basis of functions QA{x{L,Ti), 
determined for various temperatures for a fixed spin concentration p of the samples. 

6. The temperature dependence is determined for asymptotic values of the thermody- 
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namic quantity by substituting Al{t) and Q^er into relation 1^. 

Figures [TH show by way of example the scaling functions for correlation length ^ and 
susceptibility x, obtained for systems with spin concentrations p = 0.95 and 0.50 at various 
temperatures using the polynomial approximation in variable x (Figs. [T], [3]) and in variable 
exp {—1/x) (Figs. O HI). It can be seen from the figures that the scaling functions show a 
tendency towards a single universal curve for each spin concentration p in the entire range 
of variation of scaling variable xl- 

Figures 5 and 6 show the averaged scaling functions for the correlation length and suscep- 
tibility for various spin concentrations p, which were obtained using the polynomial approxi- 
mation in X (Fig. E]) and in exp (—1/x) (Fig. [H]). The averaged scaling functions demonstrate 
a tendency indicating the possible existence of two classes of universal critical behavior for 
the diluted Ising model with different modes of behavior for weakly {p = 0.95, 0.80) and 
strongly {p = 0.60, 0.50) disordered systems. 

Table [I contains asymptotic values of ^(T) and x(T) obtained using averaged scaling 
functions for various temperatures and spin concentrations. The errors in values of ^(T) 
and x{T) take into account statistical errors in the measured values of C,l{T) and xl{T), as 
well as approximation errors. 



TABLE I: Asymptotic values of correlation length and susceptibility obtained using scaling func- 
tions with a polynomial dependence on x (pol) and exp(— (exp). 
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IV. CALCULATION OF CRITICAL CHARACTERISTICS 

Asymptotic critical exponent of a thermodynamic quantity A{t) is described by the 
expression 

'^ = -li^nT4T^' A{t) = AM-^ (13) 
r^o ln|r| 

where and A- are the critical amplitudes above and below the critical temperature, 
respectively. A power law of the type ( |T3i) is accurate only in the limit r — 0. To calculate 



critical exponents in the intermediate nonasymptotic regime, we must introduce additional 
correcting terms to power law (fT^ . In accordance with the Wegner expansion [6| we have 

A{t) = (Ao + AiT^' + A2r2-- + . . . ) r-^ (r > 0), (14) 

where Ai are nonuniversal amplitudes and uj is the critical exponent of the correction to 
scaling. Here, in the calculation of the characteristics of disordered systems, we use the first 
correction to the asymptotic behavior for the correlation length and susceptibility: 

e(r) = T-'^ (4 + Ar') , e = u;u, (15) 

X(r) = T-^ {A^ + A^y) , (16) 

and calculated critical exponents z/, 7 and 6, as well as the critical temperatures using 
the least squares method for the best approximation of the data presented in Table [T] by 
expressions f|T5l) and (|T6l) . Table HTl contains the values of critical parameters obtained for 
various spin concentrations p using the initial data corresponding to various approximations 
for scaling functions, as well as their values averaged over approximations. It can be seen 
that the critical exponents form two groups with close values to within experimental error. 
The first group corresponds to p = 0.95 and 0.80, i.e., to weakly disordered systems with 
spin concentrations p, larger than the impurity percolation threshold pimp {pimp = 0.69 for 
cubic systems), while the second group with p = 0.60 and 0.50, corresponds to strongly 
disordered systems with p^ < p < Pimp, where Pc is the spin percolation threshold {pc = 0.31 
for cubic systems) , ; in the latter case, two mutually penetrating (spin and impurity) clusters 
exist in the system. Fractal effects of these two penetrating clusters may be responsible for 
the change in the type of critical behavior for strongly disordered systems. We can consider 
that the averaged values of critical exponents u = 0.693(5), 7 = 1.342(7) and 9 = 0.157(92) 
for weakly disordered systems and u = 0.731(11), 7 = 1.422(12) and 6 = 0.203(106) for 
strongly disordered systems are the final results of our investigations. It should be noted 
that the values of the critical exponents obtained for weakly disordered systems correlate 
with the values of = 0.678(10), 7 = 1.330(17) and = 0.170(71) (to = 0.25(10)), obtained 
in [3] by the renormalizations group methods in the six-loop approximation, which are valid 
only for systems with low concentrations of impurities. 

The above values of critical exponents u and 7 are also in good agreement with the 
available results of experiments with diluted Ising-like magnets (Table UTTl) . 



TABLE II: Values of critical parameters for two types of approximations (pol) and (exp) and their 
averaged (aver) values for systems with various spin concentrations p 







p=0.95 pol 
exp 
aver 


0.6883(15) 1.3339(25) 0.141(52) 0.152(50) 4.26264(4) 4.26269(3) 
0.6935(26) 1.3430(33) 0.113(64) 0.142(54) 4.26265(5) 4.26270(3) 
0.6909(33) 1.3385(54) 0.137(56) 4.26267(4) 


p=0.80 pol 
exp 

aver 


0.6960(29) 1.3473(30) 0.180(107) 0.193(74) 3.49937(21) 3.49954(14) 
0.6947(28) 1.3421(30) 0.147(94) 0.192(71) 3.49940(21) 3.49961(14) 
0.6956(29) 1.3447(40) 0.178(87) 3.49948(18) 


p=0.60 pol 
exp 
aver 


0.7272(37) 1.4253(34) 0.221(147) 0.201(63) 2.42409(11) 2.42404(6) 
0.7233(24) 1.4054(43) 0.184(92) 0.192(109) 2.42414(8) 2.42423(7) 
0.7253(36) 1.4154(107) 0.199(103) 2.42413(9) 


p=0.50 pol 
exp 
aver 


0.7372(25) 1.4299(26) 0.164(159) 0.195(74) 1.84503(7) 1.84512(3) 
0.7368(26) 1.4266(30) 0.242(96) 0.226(66) 1.84503(7) 1.84519(3) 
0.7370(33) 1.4283(33) 0.207(100) 1.84509(6) 



TABLE III: Values of critical exponents and 7 experimentally measured by different authors in 
materials corresponding to the disordered Ising model 
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0.75(5) 


1.57(16) 



TABLE IV: Values of critical exponents v and 7 obtained by different authors using the Monte 
Carlo method 



Authors 
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Wang et al., 1990, [13] 


300 
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1.36(4) 


Heuer, 1993, [14] 
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0.95 


0.64(2) 
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Wiseman et al.,1998, [15] 
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1.357(8) 
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1.508(28) 


Ballesteros et al., 1998, [16] 


128 


0.4 < p < 0.9 0.6837(53) 


1.342(10) 0.253(43) 


Calabrese et al, 2003, [17] 


256 


0.8 


0.683(3) 


1.336(8) 0.581(85) 


Berche et al, 2005, [18] 


96 


0.85 


0.662(2) 


1.314(4) 


Murtazaev et al., 2004, [19] 


60 


0.95 


0.646(2) 


1.262(2) 






0.9 


0.664(2) 
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Table [IV] contains the latest results of Monte Carlo simulation of the critical behavior of 
the diluted Ising model obtained by various authors. Each work cited in the table has its 
merits due to the application of various techniques for processing the results of simulation, 
as well as disadvantages associated with the small size of experimental lattices, which does 
not ensure reliable asymptotic values of the quantities being measured, or with insufficient 
statistics of averaging over various impurity configurations for obtaining reliable results, 
or with disregard of the effect of nonasymptotic corrections to scaling in the calculation 
of critical exponents (the inclusion of these corrections is especially important for samples 
with spin concentrations p = 0.95 and 0.90 and strongly disordered systems). The results 
obtained in p^, 119|] can be treated as supporting our conclusions; these ideas were formulated 



in our earlier publications [20] on computer simulation of critical dynamics of the disordered 
Ising model. In fact, in spite of the attractiveness of the idea about a single universal 
critical behavior with the asymptotic values of critical exponents independent of the spin 



concentration, which was supported by the authors of the results obtained in [16| did 
not make it possible to adequately explain the results obtained for samples with p = 0.90 
using the universal critical exponent of scaling correction uj = 0.37(6) for all systems. At 
the same time, the nonasymptotic values of critical exponents obtained in 16| demonstrated 
explicit dependence on p and, assuming that u is not unique, led to two sets of asymptotic 
critical exponents for weakly and strongly disordered systems. The results of the remaining 
studies carried out on samples with a single spin concentration coincide as a rule with our 
results to within experimental error, although some mismatching associated in all probability 
with the above-mentioned drawbacks also takes place. 



V. CONCLUSIONS 



The results of our investigations lead to the following conclusions. 

(i) Scaling functions and values of critical exponents for the correlation length and suscep- 

tibility demonstrate the existence of two classes of universal critical behavior for the 
diluted Ising model with various characteristics for weakly and strongly disordered 
systems. 

(ii) The values of critical exponents obtained for weakly disordered systems are in good 

agreement with the results of the field-theoretical description to within statistical 
errors of simulation and the numerical approximations used. 

(iii) The results of simulation match the results of experimental studies of the critical 
behavior of diluted Ising-like magnets. 
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